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1.  Introduction 


One  of  the  major  challenges  facing  researchers  studying  nanoscale  systems  is  iden¬ 
tifying  the  dynamic  mechanisms  controlling  material  properties  and  performance. 
This  is  especially  apparent  in  nanocrystalline  metals  where  mechanical  properties 
are  largely  governed  by  atomic  scale  defects  and  interfaces  between  crystalline  re¬ 
gions.  Molecular  dynamics  simulations  provide  researchers  with  a  tool  to  accu¬ 
rately  simulate  such  small  systems  with  atomic  resolution,  but  these  simulations 
lack  an  efficient  means  of  producing  the  complex  structures  associated  with  nano¬ 
grained  systems.  These  limitations  mean  that  researchers  must  generate  nanocrys¬ 
talline  structures  manually  with  an  external  algorithm. 

In  this  technical  note,  an  algorithm  (nanocry staljbuilder.py)  implemented  in  Python 
is  used  to  efficiently  generate  highly  tailored  nanocrystalline  microstructures  based 
on  user  input.  The  Python  script  is  attached  as  the  Appendix  and  a  description  of  its 
capabilities  and  execution  is  described  within  this  document. 

2.  Background 

2.1  Voronoi  Tessellation 

The  nanocrystalline  builder  algorithm  is  based  on  a  Voronoi  tessellation  method  that 
implements  packing  rules  to  create  optimal  grain  morphologies.  Voronoi  tessella¬ 
tions,  also  called  Voronoi  diagrams,  are  powerful  space-filling  geometric  structures 
that  have  many  applications  in  modern  science.  In  this  work,  center  of  mass  Voronoi 
tessellations  are  used  to  construct  3-dimensional  polyhedra,  which  define  the  limits 
of  individual  grains  in  a  polycrystal.  In  center  of  mass  Voronoi  tessellations,  a  set  of 
points,  p,  are  defined  as  local  sources  of  influence  (centers  of  mass)  for  a  spherical 
field  q.  Individual  points  q;  are  assigned  to  centers  of  influence  p,  based  on  their 
proximity.  Voronoi  regions  are  then  defined  as  the  collection  of  all  points  that  are 
assigned  to  a  single  region  of  influence.  Figure  1  showcases  a  simple  2-dimensional 
case  with  25  centers  of  mass,  where  the  Voronoi  regions  are  bounded  by  the  center 
lines  separating  Voronoi  centers. 

Traditionally,  Voronoi  tessellation  algorithms  are  initialized  with  either  randomly 
placed  center  of  masses  or  centers  packed  using  a  regular  grid;  however,  these  are 
not  ideal  for  constructing  nanocrystals.  The  disadvantage  with  implementing  ran¬ 
domly  dispersed  Voronoi  tessellation  algorithms  for  nanocrystalline  construction  is 
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that  if  the  Voronoi  (grain)  centers  are  not  properly  packed,  the  resulting  structure 
will  likely  be  unrealistic.  Purely  random  placement  of  centers  can  result  in  acute 
grains  and  needle-like  grain  structures  that  are  rarely,  if  ever,  observed  in  physical 
systems.  Similarly,  grid  packing  can  result  in  highly  cubic  grains  or  idealized  grains 
that  are  equally  unlikely.  Therefore,  in  this  work,  a  spherically  optimized  packing 
routine  was  implemented  to  overcome  the  hurdles  presented  by  these  traditional 
packing  routines  to  create  more  realistic  nanocrystalline  starting  structures. 


Fig.  1  A  Voronoi  diagram  showing  the  tessellation  for  25  randomly  generated  centers  of  mass. 
The  centers  are  blue  circles.  The  blue  lines  are  equal  distance  to  the  2  closest  center  points. 
The  red  triangles  are  vertices  points  at  equal  distance  to  3  center  points. 


2.2  Nearest-Neighbor  Analysis 

One  obstacle  to  using  Voronoi  tessellations  for  creating  atomic  nanocrystalline 
structures  is  the  computationally  inefficient  process  of  testing  for  nearest  neigh¬ 
bors  when  assigning  grain  centers  and  populating/trimming  atoms.  The  brute  force 
method  loops  over  each  point  and  repeatedly  calculates  the  distances  to  all  other 
points  (incorporating  periodic  boundaries)  to  define  nearest  neighbors.  While  this 
method  is  easy  to  implement,  such  brute  force  methods  waste  computational  re¬ 
sources  by  running  calculations  on  nonlocal  points  that  could  be  ignored.  It  is 
possible  to  include  rules  for  ignoring  points  far  away,  which  can  help  this  brute 
force  approach  to  be  slightly  more  efficient.  However,  since  nanocrystal_builder.py 
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generates  potentially  millions  of  data  points  (atoms),  a  more  efficient  method  for 
proximity  testing  was  necessary.  Typically,  high  (computational)  efficiency  nearest- 
neighbor  analysis  is  performed  by  breaking  the  data  set  down  into  a  tiered  data 
structure  that  spatially  sorts  the  atom  coordinates  into  bins.  This  allows  these  high- 
efficient  algorithms  to  limit  the  number  of  points  queried  to  just  those  belonging  to 
the  same  tier  (or  neighboring  tiers)  as  the  point  of  interest.  An  example  of  such  a 
data  structure,  the  KDTree  from  SciPy  stack  in  Python,  is  leveraged  in  this  work  to 
improve  the  computational  efficiency  of  the  nearest-neighbor  testing  routines. 

3.  Methodology 

The  following  sections  provide  a  detailed  description  of  how  nanocry staljbuilder.py 
works.  The  source  code  has  been  broken  apart  into  several  subfunctions  to  help  im¬ 
prove  readability  and  to  establish  a  modular  flow  control.  Each  function  in  the  code 
has  been  designed  to  perform  a  specific  set  of  tasks  that  will  be  described  separately. 

In  general,  there  are  2  ways  to  construct  nanocrystals  with  nanocry stal_builder.py . 
The  first  method  uses  an  optimized  packing  algorithm  to  generate  a  list  of  grain 
centers  that  are  populated  with  seeds — spherical  groups  of  atoms  extracted  from  a 
reference  file.  This  method  uses  a  single  reference  file,  and  is  therefore  referred 
to  as  Single  Mode.  The  second  method,  referred  to  as  Config  Mode,  uses  data 
from  a  configuration  file  to  generate  nanocrystalline  structures  thus  providing  addi¬ 
tional  control  over  the  structure.  Configuration  files  contain  information  pertaining 
to  the  simulation  cell  size,  average  grain  size,  as  well  as  the  position,  orientation 
and  reference  data  for  each  grain.  Configuration  files  are  automatically  generated 
by  nanocry stal_builder.py  as  a  standard  output  to  ensure  the  resulting  structure  is 
easily  reproducible.  Users  can  modify  configuration  files  to  obtain  a  high  degree 
of  control  over  the  resulting  structure,  allowing  for  the  construction  of  networked 
grains  and  rescaled  grain  sizes. 
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3.1  Required  Software 


The  nanocrystaljbuilder.py  algorithm  detailed  in  this  document  is  written  in  Python 
and  makes  use  of  the  SciPy  stack.  Development  and  testing  were  performed  using 
Python  2.6.6  and  the  SciPy  0.15.0  with  its  associated  libraries.  These  version  num¬ 
bers  should  be  considered  the  minimum  required  for  use.  The  standard  output  for 
the  atomic  data  generated  by  nanocrystaljbuilder.py  is  formatted  to  match  the  data 
file  format  associated  with  the  molecular  dynamics  code  LAMMPS.1 2  Therefore, 
structures  can  be  read  into  LAMMPS  using  the  read_data  command  (i.e., 
http://lammps.sandia.gov/doc/read_data.html). 

3.2  Single  Mode  (Generating  from  a  LAMMPS  Dump  File) 

This  section  discusses  the  methods  and  code  used  to  generate  a  nanocrystalline 
structure  with  a  single  reference  file  for  seed  extraction.  Some  of  the  code  segments 
detailed  here  are  also  used  in  Config  Mode  as  discussed  in  Section  3.3.  It  is  recom¬ 
mended  that  the  reference  structures  used  for  Single  Mode  be  of  single  crystalline 
or  bicrystalline  nature  as  more  complex  structures  may  result  in  unrealistic  struc¬ 
tures.  Additional  care  must  be  taken  when  using  a  bicrystalline  reference  structure 
as  the  existence  of  a  grain  boundary  within  the  seed  effectively  reduces  the  grain 
size  and  may  create  potentially  unrealistic  nanocrystalline  structures. 

3.2.1  Wrapping  Periodic  Bounds 

The  nanocrystaljbuilder.py  script  accommodates  user-defined  periodic  boundaries. 
This  task  is  performed  by  testing  if  a  point,  which  can  be  an  atom  or  a  grain  center, 
is  positioned  beyond  the  upper  and  lower  bounds  of  the  simulation  cell  for  each 
dimension.  If  one  of  the  boundaries  is  exceeded,  the  point  is  conditionally  moved 
by  adding  (or  subtracting)  a  periodic  length  if  the  point  lies  outside  of  the  lower 
(or  upper)  bound.  Performing  this  process  in  all  3  axes  results  in  all  points  being 
contained  within  the  simulation  cell. 

3.2.2  Optimized  Packing  of  Grain  Centers 

The  nanocrystaljbuilder.py  script  uses  Voronoi  tessellations  to  form  3 -dimensional, 
grain-like  structures.  In  this  routine  an  initial  center  is  placed  at  the  corner  of  the 
simulation  cell  lying  on  the  origin.  A  cloud  of  allowed  points  is  then  defined  at  a 
minimum  and  maximum  radius  around  this  center;  these  radii  are  defined  to  be  75% 
and  125%  of  the  requested  grain  size,  respectively.  The  next  grain  center  is  chosen 
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at  random  from  these  points  and  a  new  cloud  of  allowed  points  is  defined  such  that 
none  of  the  new  points  fall  within  the  minimum  radius  of  any  existing  grain  center. 
This  method  is  repeated  until  an  iteration  is  reached  in  which  the  number  of  allowed 
points  becomes  zero,  thereby  optimally  filling  the  requested  simulation  cell  size.  To 
account  for  periodic  boundary  conditions,  “ghost”  grain  centers  are  placed  in  the 
area  outside  the  simulation  cell  which  represent  the  reflections,  across  the  periodic 
bounds,  of  previously  defined  centers.  This  ensures  that  centers  near  the  simulation 
edge  are  not  within  the  minimum  radius  of  each  other  when  wrapped  around  the 
periodic  boundaries. 

3.2.3  Generating  and  Scaling  Reference  Files 

The  nanocry stal_builder.py  script  uses  reference  files  constructed  in  the  LAMMPS 
dump  file  format.  In  the  event  that  the  reference  file  provided  by  the  user  is  smaller 
than  the  desired  grain  size,  the  nanocry staljbuilder.py  code  enters  into  a  dynamic 
rescaling  function.  This  function  checks  the  dimensions  of  the  reference  file  and 
replicates  the  reference  atoms  across  the  periodic  edge.  This  is  repeated  until  all  of 
the  reference  file  dimensions  are  at  least  150%  of  the  requested  grain  size.  In  Single 
Mode  this  is  done  once  and  the  newly  replicated  reference  is  held  in  memory  until 
the  code  exits.  In  Config  Mode  the  reference  structure  for  each  grain  is  individually 
rescaled  and  dropped  from  memory  at  the  end  of  each  iteration  of  the  code. 

3.2.4  Populating  Atoms 

The  nanocry  staljbuilder.py  script  then  populates  each  Voronoi  grain  with  atoms 
from  the  rescaled  reference  file.  Figure  2  highlights  the  steps  that  are  used.  Images 
generated  using  OVITO.3  The  first  step  in  this  process  is  to  isolate  a  sphere  of  atoms 
from  the  reference  file.  This  is  done  by  searching  for  all  atoms  within  a  cutoff  radius 
of  the  reference  structure  center  point  using  a  KDTree  query.  By  default  the  sphere 
of  atoms  is  rotated  randomly  using  the  axis-angle  method.  However,  if  a  texturing 
axis  has  been  specified,  the  rotation  axis  is  fixed  to  the  texturing  axis.  The  rotated 
sphere  of  atoms  is  then  assigned  to  a  grain  center  by  placing  the  center  point  of  the 
sphere  on  the  grain  center.  This  process  is  repeated  until  all  of  the  grain  centers  are 
populated. 

At  this  point,  the  structure  is  populated  with  overlapping  spheres  of  atoms.  These 
spheres  are  trimmed  according  to  the  Voronoi  tessellation  of  the  simulation  cell 
(i.e.,  the  boundaries  between  grains  are  defined  by  the  midlines  between  grain  cen- 
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ters  [see  Fig.  1]).  To  determine  if  an  atom  lies  within  the  boundaries  of  its  parent 
grain,  KDTree  queries  are  used  to  test  if  the  nearest  grain  center  to  the  atom  is  the 
grain  center  it  was  assigned  (its  parent).  If  the  atom  is  closest  to  a  center  that  is  not 
its  parent  center,  it  is  deleted.  To  account  for  the  periodic  edges  “ghost”  centers  are 
placed  around  the  defined  simulation  cell.  After  overlapping  atoms  are  removed, 
the  structure  is  wrapped  across  the  periodic  bounds  to  ensure  that  all  atoms  are 
contained  within  the  user-defined  simulation  cell  as  shown  in  Fig.  2. 


Fig.  2  Images  showing  steps  taken  when  generating  nanocrystals  (left  to  right):  populating 
cell  with  grain  centers,  sphere  of  atoms  with  defined  crystal  structure  centered  at  each  grain 
center,  identifying  atoms  belonging  to  each  grain  center,  removing  overlapping  atoms,  and 
wrapping  atoms  through  periodic  boundaries  back  into  the  simulation  cell 


3.3  Config  Mode:  Generating  from  a  Configuration  File 

The  nanocry st al_builder.py  script  can  also  operate  in  Config  Mode  to  generate 
nanocrystalline  structures.  In  this  mode,  the  nanocrystalline  structure  is  generated 
using  configuration  files  that  contain  the  simulation  cell  dimensions,  average  grain 
size,  grain  positions,  grain  orientations,  and  reference  data.  Config  Mode  enables 
easy  reproduction  of  previously  built  nanostructures  and  provides  capabilities  of 
rapidly  modifying  existing  nanostructures.  For  example,  Config  Mode  allows  the 
user  to  automatically  replace  a  percentage  of  the  grain  reference  files,  reorient 
grains  to  add  texturing,  and  to  rescale  the  entire  system.  These  options  are  enabled 
using  a  series  of  questions  posed  to  the  user  at  the  command  prompt. 

3.3.1  Reading  Grain  Information 

The  nanocry staljbuilder.py  script  reads  grain  center  information  from  a  configura¬ 
tion  file,  prompts  the  user  for  modifications,  and  then  builds  the  structure  according. 
An  example  configuration  file  is  shown  in  Fig.  3.  The  first  line  of  the  configuration 
file  is  for  user  comments  about  the  structure,  which  is  ignored  by  the  code.  The  sec¬ 
ond  line  consists  of  7  entries  pertaining  to  the  average  grain  size  and  simulation  cell 
limits,  specifically  grain  size  (A),  A min,  A j!f (:r ,  Ymin,  Ymax,  Y  !f) rr/ ,  and  Zmax.  The 
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remaining  lines  contain  information  regarding  individual  grains  in  the  following 
order:  grain  number,  center  x,  y,  z  coordinate,  rotation-axis  x,  y,  z  index,  rotation 
angle,  reference  file.  To  improve  efficiency,  the  grain  data  are  stored  as  Python  dic¬ 
tionary  objects. 


#data  for  centroids 

50.000000  0.000000  100.000000  0.000000  100.000000  0.000000  100.000000 
0  25.000000  25.000000  25.000000  0.577350  0.577350  0.577350  1.019442  out.crystcufix.txt 

1  15.653842  14.213959  73.605896  0.577350  0.577350  0.577350  2.977697  out.twinlc.txt 

2  10.845016  80.767845  21.027548  0.577350  0.577350  0.577350  0.997454  OUt.tWinlC.txt 

3  35.462112  77.228439  75.485779  0.577350  0.577350  0.577350  1.959592  out.crystcufix.txt 

4  72.080017  5.680593  99.983498  0.577350  0.577350  0.577350  2.679047  out.twinlc.txt 

5  48.061413  66.311081  38.980812  0.577350  0.577350  0.577350  1.239753  out.crystcufix.txt 

6  77.514809  2.187867  55.050816  0.577350  0.577350  0.577350  2.508741  out.twinlc.txt 

7  77.514809  45.743797  67.839989  0.577350  0.577350  0.577350  1.082618  out.crystcufix.txt 

8  39.543736  37.302536  91.937027  0.577350  0.577350  0.577350  1.483817  out.crystcufix.txt 

9  77.514809  45.743797  15.866305  0.577350  0.577350  0.577350  1.644024  out.twinlc.txt 

10  90.813972  77.228439  83.444176  0.577350  0.577350  0.577350  1.879738  out.twinlc.txt 

11  11.178849  49.204056  51.654900  0.577350  0.577350  0.577350  1.983781  out.crystcufix.txt 


Fig.  3  An  example  of  the  contents  of  a  config  file:  the  first  line  is  a  comment  line  and  is  ignored 
by  the  algorithm,  the  second  line  contains  information  about  the  grain  size  and  the  simulation 
cell  size  (listed  as  grain  size  (A),  Xmin,  Xmax,  Ymin,  Ymax,  Zmin,  and  Zmax),  and  the  remaining 
lines  contain  the  grain  center  ID,  atomic  Cartesian  coordinates,  rotation  axis,  rotation  angle 
(radians),  and  reference  file 


3.3.2  Repopulating  Atoms 

In  Config  Mode,  the  nanocry stal_builder.py  script  populates  atoms  using  the  infor¬ 
mation  contained  within  the  configuration  file  as  well  as  any  modifications  specified 
by  the  user  in  the  command  prompt.  First,  the  nanocry staljbuilder.py  script  assesses 
any  user-specified  modifications.  This  starts  by  multiplying  all  lengths  and  center 
coordinates  by  the  user-specified  rescale  factor.  Once  the  data  have  been  rescaled 
the  code  determines  how  many,  if  any,  of  the  grains  need  to  have  their  reference 
file  replaced  to  satisfy  the  user-specified  secondary  reference  composition.  Finally, 
if  the  user  has  specified  a  texturing  direction,  the  code  will  replace  the  rotation  axis 
indices  for  all  entries  with  those  requested  by  the  user. 

Once  the  grain  center  data  have  been  read  and  processed,  the  atom  population  sub¬ 
function  commences  in  the  same  manner  as  Single  Mode.  However,  as  explained 
earlier,  while  in  Config  Mode  the  code  reads  and  rescales  the  reference  file  at  each 
step  of  the  algorithm  to  accommodate  the  use  of  multiple  references.  This  additional 
processing  means  that  Config  Mode  runs  slower  than  Single  Mode,  which  performs 
the  reference  rescaling  only  once. 
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4.  Usage 


The  following  sections  describe  how  to  invoke  and  interact  with  the  nanocrys- 
tal_builder.py  script.  These  sections  act  as  a  general  user  guide  to  the  nanocry  s- 
taljbuilder.py  script. 

4.1  Invoking  the  User  Interface 

The  nanocry stal_builder.py  script  is  written  as  a  Python  run  script  and  can  therefore 
be  invoked  through  either  a  bash  terminal  or  a  Python  interactive  environment.  Prior 
to  invoking  the  script,  the  user  must  change  the  present  working  directory  to  the 
directory  containing  the  nanocry staljbuilder.py  script.  Once  this  is  done  the  script 
is  invoked  by  entering  python  nanocrystal_builder . py  in  terminal  or 
import  nanocry st al_builder  .py  in  the  Python  interactive  environment. 
Once  invoked,  the  algorithm  begins  querying  the  user  for  necessary  information. 

The  nanocry  staljbuilder.py  script  relies  on  user  interaction  to  establish  key  ini¬ 
tialization  parameters  within  the  algorithm.  This  interaction  takes  place  through  a 
series  of  simple  text  queries  to  the  user  within  the  run  environment.  The  first  query 
“single  or  config?”  sets  the  variable  environment  and  determines  what  set  of  queries 
the  user  is  subsequently  prompted  to  answer.  This  prompt  only  accepts  user  feed¬ 
back  of  “single”  or  “config”  (no  quotes).  Any  other  input  will  result  in  an  error 
prompting  the  user  to  enter  a  valid  response. 

Single  Mode.  If  the  user  responds  to  the  first  query  with  “single”  the  following 
prompts  will  be  provided: 

•  “Input  reference  file  path”: 

The  answer  to  this  query  tells  the  algorithm  where  the  crystal  reference  file  is 
located.  The  script  expects  the  user  to  return  a  valid  file  path. 

•  “Input  desired  grain  size  in  angstroms”: 

User  should  respond  with  the  desired  average  grain  size,  in  angstroms,  for 
their  structure.  Response  should  be  an  integer  or  float  value. 

•  “Input  desired  box  dimensions  in  angstroms. 

Use  the  following  convention  (xlo  xhi  ylo  yhi  zlo  zhi)”: 
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The  user  should  provide  the  desired  simulation  cell  bounds.  Individual  values 
should  be  in  integer  or  float  format  with  a  single  space  in  between  values. 
Note  that  no  parenthesis  should  be  present  in  the  response. 

•  “Input  crystalline  texturing  direction 

Reply  “none”  (quotations  excluded)  if  no  texturing  desired 

Use  following  convention  “index_x,  index_y,  index_z”: 

If  the  user  desires  their  structure  to  possess  crystalline  texturing,  they  should 
respond  to  this  query  with  the  desired  texturing  direction  using  the  Miller  in¬ 
dex  convention.  Individual  indices  should  be  integers  with  commas  separat¬ 
ing  indices.  If  texturing  is  not  desired,  the  user  should  respond  with  “none” 
(no  quotes). 

This  mode  outputs  a  configuration  file,  which  can  be  used  to  recreate  and  manipu¬ 
late  the  structure  through  Config  Mode,  and  a  LAMMPS-style  data  file  that  contains 
the  atom  types  and  positions. 

Config  Mode.  If  the  user  responds  to  the  first  query  with  “config”  (no  quotes)  the 
following  prompts  will  be  provided: 

•  “Input  path  of  configuration  file”: 

The  user  should  provide  a  valid  path  to  a  file  containing  the  desired  reference 
data.  Information  on  how  this  file  is  created  and  structured  is  provided  in  the 
methods  section. 

•  “Input  a  scaling  factor  (1  if  no  change  desired,  must  be  greater  than  0)”: 

This  query  is  used  to  linearly  rescale  the  reference  file.  Values  larger  than  1 
will  result  in  an  increase  in  size  while  values  less  than  one  will  result  in  a 
decrease  in  size.  A  value  of  1  will  result  in  no  change  in  size.  User  should 
provide  a  response  in  either  integer  or  float  format. 

•  “Percent  of  grains  to  be  reassigned  (decimal  from  0-1)”: 

User  will  specify  what  percentage  of  the  grains  in  their  structure  will  have 
their  reference  file  replaced.  Input  should  be  a  float  between  0  and  1  corre¬ 
sponding  to  the  decimal  form  of  the  desired  percentage. 
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•  “Input  path  of  new  reference  file”: 

Conditional  query,  only  requested  if  previous  is  non-zero.  User  must  specify 
a  valid  path  to  the  desired  secondary  reference  file. 

•  “Input  crystalline  texturing  direction 

Reply  “none”  (quotations  excluded)  if  no  texturing  desired 

Use  following  convention  “index_x,  index_y,  index_z”: 

If  the  user  desires  their  structure  to  possess  crystalline  texturing,  they  should 
respond  to  this  query  with  the  desired  texturing  direction  using  the  Miller  in¬ 
dex  convention.  Individual  indices  should  be  integers  with  commas  separat¬ 
ing  indices.  If  texturing  is  not  desired,  the  user  should  respond  with  “none” 
(no  quotes). 

Config  Mode  will  produce  a  new  configuration  file  containing  all  information  perti¬ 
nent  to  the  newly  created  structure  along  with  a  LAMMPS- style  data  file  containing 
atom  information. 

Final  Input  and  Runtime.  Once  the  user  has  defined  the  algorithm-critical  vari¬ 
ables  the  final  prompt  asks  the  user  for  an  output  name.  This  output  should  be  the 
general  name  as  the  appropriate  file  extensions  and  flags  are  then  automatically 
added.  At  this  point  the  algorithm  executes  and  generates  the  nanocrystalline  struc¬ 
ture.  At  various  stages  during  runtime  the  algorithm  prints  status  messages  to  the 
screen  to  verify  its  progress.  Since  default  Python  runtimes  do  not  take  advantage  of 
parallel  computing,  large  structures  can  take  some  time  to  generate  and  may  require 
significant  memory.  In  testing,  structures  containing  up  to  5  million  atoms  could  be 
generated  with  fewer  than  4  GB  of  RAM  (random  access  memory). 

4.2  Output 

Two  files  are  output  by  the  nanocry staljbuilder.py  script  upon  completion.  The  first 
file  is  the  structure  configuration  file.  As  previously  detailed,  this  file  contains  in¬ 
formation  that  can  be  used  to  reproduce  and  modify  the  structure  while  retaining 
the  overall  grain  morphology.  This  file  will  be  named  filename. config  where  file¬ 
name  is  the  user-specified  base  file  name.  The  second  file,  named  filename. data, 
is  a  LAMMPS  data  file  containing  the  ID,  type,  and  position  of  all  of  the  atoms 
in  the  structure.  By  default,  the  atom  types  will  be  based  on  which  grain  the  atom 
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belongs  to.  The  preamble  in  this  file  is  standard  to  the  LAMMPS  data  file  format 
and  contains  information  about  the  simulation  cell  dimensions. 

5.  Examples 

In  this  section,  both  Single  Mode  and  Config  Mode  for  generating  nanoscale  struc¬ 
tures  using  nanocrystaljbuilder.py  are  demonstrated.  Example  command  prompt 
options  for  these  modes  are  shown  by  the  screen  shot  in  Fig.  4.  First,  a  10-nm  grain 
structure  is  created  in  a  15-  x  15-  x  15-nm  simulation  cell.  Here,  each  grain  con¬ 
tains  an  ideal  face-centered  cubic  structure  as  described  by  the  dump.Cu  FAMMPS 
file  (that  is  located  in  the  working  directory).  No  texture  direction  is  specified,  thus 
each  grain  has  a  random  orientation.  Assuming  100%  packing  efficiency,  the  algo¬ 
rithm  estimates  that  six  10-nm  grains  need  to  be  created  within  the  15-nm  cube; 
however,  the  optimized  spherical  packing  algorithm  is  only  able  to  successfully 
place  5  grain  centers  within  the  volume.  After  determining  the  grain  center  loca¬ 
tions,  the  nanostructure  is  built  and  the  resulting  outputs  are  a  FAMMPS  data  file 
called  data.Cu_NC  and  a  configuration  file  saved  as  Cu_NC_centroids. config. 

The  nanocrystal_builder.py  script  is  invoked  a  second  time  to  demonstrate  the  Con¬ 
fig  Mode  in  the  lower  half  of  Fig.  4.  Here,  the  Cu_NC_centroids. config  file  gen¬ 
erated  from  the  previous  example  is  used  to  preserve  the  grain  morphology  (grain 
center  locations).  However,  in  this  example  the  simulation  dimensions  are  rescaled 
1.5  times  to  create  approximately  15-nm  grains.  Using  the  command  prompt,  75% 
of  the  grains  are  randomly  chosen  to  be  reassigned  to  a  secondary  atomic  struc¬ 
ture.  The  secondary  atomic  structure,  called  dump.  Cu_twin  (located  in  the  work¬ 
ing  directory),  contains  approximately  5-nm  spaced  Cu  twins.  Through  command 
prompt  interaction,  all  grain  orientations  are  reoriented  about  the  [111]  direction 
to  observe  crystallographic  texturing.  As  previously  shown,  the  resulting  outputs 
from  the  nanocry staljbuilder.py  code  are  a  structure  and  a  configuration  file.  Both 
structures  created  in  this  section  are  shown  in  Fig.  5. 
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bash-4.1  $  python  nanocrystal_builder.py 
Do  you  want  to  generate  from  a  single  reference  or  a 
configuration  file  (respond  single  or  config): 
single 

Input  the  reference  file  path: 

./dump.Cu 

Input  desired  grain  size  in  angstroms: 

100 

Input  the  box  dimensions  in  angstroms 

Use  the  following  convention  lowerx,  upperx,  lowery, 

uppery,  lowerz,  upperz: 

0,  150,  0,  150,  0,  150 

Input  crystalline  texturing  direction 

Reply  "none"  (quotations  excluded)  if  no  texturing  desired 

Use  following  convention  index_x,  index_y,  index_z: 

none 

Input  desired  output  basename: 

Cu_NC 

asigning  grain  centers 

approximate  number  of  grains  =  6 

asigning  grain  1 

asigning  grain  2 

asigning  grain  3 

asigning  grain  4 

asigning  grain  5 

populating  atoms 

populating  grain  0 

populating  grain  1 

populating  grain  2 

populating  grain  3 

populating  grain  4 

populating  grain  5 

trimming  excess  atoms 

wrapping  periodic  bounds 

outputting  to  data  file 


bash-4. 1$  python  nanocrystal_builder.py 
Do  you  want  to  generate  from  a  single  reference  or  a 
configuration  file  (respond  single  or  config): 
config 

Input  path  of  configuration  file: 

Cu_NC_centroids. config 

Input  a  scaling  factor  (1  if  no  change  desired,  must  be 
greater  than  0): 

1.5 

Percent  of  grains  to  be  reassigned  (decimal  from  0-1): 
0.75 

Input  path  of  new  reference  file: 

,/dump.Cu_twin 

Input  crystalline  texturing  direction 

Reply  "none"  (quotations  excluded)  if  no  texturing  desired 

Use  following  convention  index_x,index_y,index_z: 

1,1,1 

Input  desired  output  basename: 
Cu_NC_1.5scale_0.75twin_111 
reading  configuration  data 
reassigning  4  reference  files 
populating  atoms 
populating  grain  0 
populating  grain  1 
populating  grain  2 
populating  grain  3 
populating  grain  4 
populating  grain  5 
trimming  excess  atoms 
wrapping  periodic  bounds 
outputting  to  data  file 


Fig.  4  User  feedback  flow  for  Single  Mode  without  texturing  (top),  and  Config  Mode  with 
texturing,  rescaling,  and  reassigning  (bottom) 
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Fig.  5  Image  demonstrating  the  outputs  from  the  examples  in  Fig.  4.  The  nanocrystalline 
structure  on  the  left  is  generated  from  Single  Mode,  the  structure  on  the  right  is  generated 
from  Config  Mode,  where  it  has  been  rescaled  to  150%  original  size  and  75%  of  the  grains 
have  been  reassigned  to  have  twins. 


6.  Extensibility 

6.1  Single  Mode 

In  the  current  version,  Single  Mode  allows  control  over  grain  size,  simulation  cell 
size,  and  crystallographic  texturing.  These  capabilities  are  well  suited  for  producing 
homogeneous  structures  with  equiaxed  grains  but  are  not  well  suited  for  producing 
grain  size  gradients  or  non-equiaxed  grain  shapes.  It  is  possible  to  add  these  capa¬ 
bilities  by  manipulating  the  centroid  placement  algorithm  to  use  a  gradient  radius 
function  or  an  elliptical  region  of  allowed  space  instead  of  the  homogeneous  spheri¬ 
cal  space.  With  the  current  distribution  this  extension  should  be  considered  as  future 
work. 

6.2  Config  Mode 

The  current  version  is  configured  to  automatically  use  2  separate  reference  files 
when  rebuilding  a  structure.  This  is  very  useful  for  adding  twinned  grains  to  an  ex¬ 
isting  granular  structure,  as  was  the  intention  of  this  project.  However,  should  the 
user  desire  to  produce  a  structure  using  additional  reference  files,  they  would  have 
to  either  manually  manipulate  the  configuration  file  or  run  Config  Mode  multiple 


Approved  for  public  release;  distribution  is  unlimited. 


13 


times.  As  both  of  these  options  are  suboptimal  for  large  structures  or  many  refer¬ 
ences,  future  versions  of  this  algorithm  can  include  capabilities  to  automatically 
process  structures  with  arbitrary  numbers  of  reference  files. 

7.  Conclusions 

The  nanocrystaljbuilder.py  script  was  developed  to  generate  nanocrystalline  struc¬ 
tures  with  flexibility  to  alter  texture,  reference  structures,  and  grain  sizes.  The  al¬ 
gorithm  uses  an  optimized  packing  routine  to  produce  realistic  grain  shapes  and 
size  distributions.  An  efficient  nearest-neighbor  matching  algorithm  is  leveraged  to 
eliminate  excess  atoms.  To  increase  functionality,  the  algorithm  is  able  to  read  in 
previously  built  structures  and  modify  them  at  user  request.  This  functionality  al¬ 
lows  the  user  to  resize,  twin,  or  alloy  a  structure  while  retaining  grain  geometry. 
In  testing,  the  algorithm  successfully  produced  a  series  of  structures  in  which  the 
grain  size,  twin  density,  and  texturing  axis  were  varied  while  the  grain  geometry 
remained  consistent. 

The  algorithm  produced  from  this  research  is  currently  being  used  to  study  the 
effect  of  nanotwins  on  the  mechanical  response  during  deformation.  This  study  is 
being  complemented  by  the  use  of  virtual  diffraction  techniques  to  directly  compare 
results  to  experiments.  Ongoing  efforts  are  being  made  to  expand  the  algorithm 
to  accommodate  more  complex  structures,  such  as  size  distributions  and  irregular 
grain  sizes. 

The  nanocrystal_builder.py  code  can  be  downloaded  by  clicking  here. 
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Intentionally  lelt  blank. 
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Appendix.  Python  Script:  nanocrystal builder.py  Function 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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#  File:  nanocrystal_builder . py 

#  Script  builds  nanocrystalline  centers  using  optimized  grain  center  placement 

#  Daniel  Foley 

#  Army  Research  Laboratory 

#  Summer  2015/2016 

#  - Initialize - 

import  os 

import  sys 

import  numpy  as  np 

from  scipy  import  spatial  as  sp 

import  random 

import  linecache 

import  math  as  mt 

reference_mode  =  raw_input ( ' Do  you  want  to  generate  from  a  single  reference  or  a  configuration  file  '  \ 

' (respond  single  or  config) :\n') 

while  ref erence_mode  !=  'single'  and  ref erence_mode  !=  'config': 

print  ("I  could  not  interpret  your  choice,  please  check  your  spelling  and  try  again.") 
ref erence_mode  =  raw_input (' single  or  config?  ') 

if  reference_mode== ' single ' : 

crref  =  raw_input (' Input  the  reference  file  path:\n') 

avg_grain_size  =  float (raw_input (' Input  desired  grain  size  in  angstroms : \n ') ) 
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box_bounds  =  raw_input (' Input  the  box  dimensions  in  angstroms\nUse  the  following  convention'  \ 

'  lowerx,  upperx,  lowery,  uppery,  lowerz,  upperz : \n ' ) 
direction  =  raw_input (' Input  crystalline  texturing  direction\nReply  "none"  (quotations  excluded) '  \ 

'  if  no  texturing  desired\nUse  following  convention  '  \ 

'index_x,  index_y,  index_z:\n') 


if  direction  !=  'none': 

direction  =  np . fromstring (direction, sep= 1 ,' ) 
direcnorm  =  np . linalg . norm (direction) 
direction  =  direction/direcnorm 


elif  ref erence_mode== ' conf ig ' : 

ref_list=raw_input (' Input  path  of  configuration  file:\n') 
rescale_f actor  =  float (raw_input (' Input  a  scaling  factor  (1  if  no  change  desired,  must  be  '  \ 

' greater  than  0 ) : \n ' ) ) 

contwin  =  float (raw_input (' Percent  of  grains  to  be  reassigned  (decimal  from  0-1) :\n')) 
if  contwin  ! =  0 . : 

twinr  =  raw_input (' Input  path  of  new  reference  file:\n') 
else : 

twinr  =  ' f oo ' 

direction  =  raw_input (' Input  crystalline  texturing  direction\nReply  "none"  (quotations  '  \ 

'excluded)  if  no  texturing  desired\nUse  following  convention  '  \ 

' index_x, index_y, index_z : \n ' ) 

if  direction  !=  'none': 

direction  =  np . fromstring (direction, sep= '  ,  ') 
direcnorm  =  np . linalg . norm (direction) 
direction  =  direction/direcnorm 
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output_name  =  raw_input (' Input  desired  output  basename : \n ' ) 

# - Define  Functions - 

def  efficient_packing (center, grain_size) : 

#defines  a  point  space  in  which  the  next  grain  center  may  be  placed 
rot  =  np . linspace ( 0 . , 360*np . pi/ 1 80 . , 45 ) . reshape ( 45 , 1 ) 
r  =  np . linspace (. 75*grain_size, 1 . 25*grain_size,  20 )  . reshape (20, 1) 
zero  =  np . zeros ( (20 , 2 ) ) 
r  =  np . concatenate (( r, zero) , axis  =  1) 
point_space  =  np . empty ( [20*45, 3] ) 

n  =  0 

to 

O  for  j  in  range (0,20) :  #this  can  probably  be  vectorized...  work  in  progress 

for  i  in  range (0,45) : 

R_y  =  np . array ( [ [np . cos (rot [ i, 0 ] ) , 0 , np . sin ( rot [i,0])], [0,1,0],  \ 

[-np .sin (rot [i,0] ) ,0,np.cos (rot [ i, 0] ) ] ] ) 
point_space [n, 0 : 3 ]  =  np . dot (r [ j , 0 : 3 ] , R_y ) 
n  =  n  +  1 

point_space_2  =  np . empty ( [20*45*45, 3] ) 
n  =  0 

for  j  in  range (0, len (point_space) ) :  #this  can  probably  be  vectorized...  work  in  progress 
for  i  in  range (0,45) : 

R_z  =  np . array ( [ [np . cos ( rot [ i, 0 ] ) , -np .sin (rot [ i, 0 ] ) , 0 ] , [np . sin ( rot [ i , 0 ] ) ,  \ 

np . cos (rot [i, 0] ) , 0] ,  [0,0,1]]) 
point_space_2 [n, 0 : 3 ]  =  np . dot (point_space [ j , 0 : 3 ] , R_z ) 
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n+1 


n  = 

point_space_2  =  point_space_2  +  center 
return  point_space_2 

def  ghost_centers (centroids,  box_limits)  : 

♦defines  grain  centers  as  seen  across  periodic  bounds  as  "ghosts" 
ghost_x  =  np . array ( [- (box_limits [0,1] -box_limits [ 0 , 0 ] ) , 0 , (box_limits [0,1] -box_limits [0,0] ) ] ) 
ghost_y  =  np . array ( [- (box_limits [1,1] -box_limits [ 1 , 0 ] ) , 0 , (box_limits [1,1] -box_limits [1,0])]) 
ghost_z  =  np . array ( [- (box_limits [2,1] -box_limits [ 2 , 0 ] ) , 0 , (box_limits [2,1] -box_limits [2,0] ) ] ) 
ghost  =  np . empty ( [len (ghost_x) *len (ghost_y) *len (ghost_z) , 3] ) 
roll  =  0 

for  i  in  range ( 0 , len (ghost_x) ) :  #this  can  probably  be  vectorized...  work  in  progress 
for  j  in  range ( 0 , len (ghost_y) ) : 

for  k  in  range (0, len (ghost_z) ) : 

ghost [roll, 0]  =  ghost_x[i] 

ghost [roll, 1]  =  ghost_y[j] 
ghost [roll , 2 ]  =  ghost_z[k] 
roll  =  roll  t  1 

if  len  (np . shape (centroids ) )  ==  1: 

grain_centers_l  =  centroids . reshape ( 1 , 3 ) 
grain_centers_2  =  centroids . reshape ( 1 , 3 ) 
else : 

grain_centers_l  =  centroids 
grain_centers_2  =  centroids 
for  i  in  range ( 0 , len (ghost )) : 

if  ghost [i,0]  ==  0  and  ghost [i,l]  ==  0  and  ghost[i,2]  ==  0: 
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continue 
else : 

hold_center  =  np . add (grain_centers_l [ : , 0 : 3] ,  ghost [i, :]) 
grain_centers_2  =  np . concatenate ( (grain_centers_2 , hold_center ) , axis=0 ) 
return  grain_centers_2 

def  grain_centroid (amorphous_edge, grain_size) : 

#assigns  grain  centers  within  an  allowed  space  defined  by  the  positions  of  existing  grain  centers 
box_vol  =  (amorphous_edge [0,1] -amorphous_edge [0,0])* (amorphous_edge [1,1] -amorphous_edge [1,0])  \ 

* ( amorphous_edge [2,1] -amorphous_edge [2,0]) 
grain_volume  =  ( ( 4 . /3 . ) * (grain_size  /  2)**3)*np.pi 
N  =  int (box_vol  /  grain_volume) 

print  'approximate  number  of  grains  =  %i '  %  (N) 

centroids  =  np . array ([[ (grain_size  /  2 . ) , (grain_size  /  2 . ) , (grain_size  /  2.)]]) 

allowed  =  efficient_packing (centroids [0, 0 : 3] , grain_size) 

allowed  =  transpose_periodic_bounds (allowed, amorphous_edge) 

ghost  =  ghost_centers (centroids [0 ,:], amorphous_edge) 

for  i  in  range ( 1 , N) : 

print  'asigning  grain  %i'  %  (i) 

if  len (allowed)  ==  0: 

print  'out  of  space' 
break 

choose  =  random. randint (0, len (allowed) -1) 

centroids  =  np . concatenate ( (centroids , allowed [ choose, 0:3]. reshape (1,3) ) ) 

allowed  =  np . concatenate ( (allowed, ef f icient_packing (centroids [ i, 0 : 3 ] , grain_size) ) ) 

cent tree  =  sp . cKDTree (ghost_centers ( centroids [ 0 : i  +  1 ,  : ] , amorphous_edge) ) 
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dist,  index  =  centtree . query (allowed) 

deny  =  np . where (dist  <  . 75*grain_size, False, True) 

if  np. any (deny)  ==  'False': 

print  'out  of  space,  %i  grains  assigned'  %  (i) 
break 

allowed  =  allowed [deny, : ] 

allowed  =  transpose_periodic_bounds (allowed, amorphous_edge) 
return  centroids 


def  seed_rotations (refdat, grain_center, grabradius, grabedges, axis, angle) : 

#function  extracts  and  rotates  crystalline  seeds  (updated  in  version  2) 
holdxl  =  int  (grabedges [ 0 , 0 ] +grabradius ) 
holdx2  =  int (grabedges [ 0 , 1 ] -grabradius ) 
if  holdxl  ==  holdx2 : 

holdxl  =  holdxl  -  2 
holdx2  =  holdx2  +  2 
if  holdxl  >  holdx2 : 

holdxl  =  int (grabedges [ 0 , 1 ] -grabradius ) 
holdx2  =  int (grabedges [ 0 , 0 ] tgrabradius ) 
holdyl  =  int (grabedges [ 1 , 0 ] igrabradius ) 
holdy2  =  int (grabedges [ 1 , 1 ] -grabradius ) 
if  holdyl  ==  holdy2 : 

holdyl  =  holdyl  -  2 
holdy2  =  holdy2  +  2 
if  holdyl  >  holdy2: 
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holdyl  =  int (grabedges [ 1 , 1 ] -grabradius ) 
holdy2  =  int (grabedges [ 1 , 0 ] +grabradius ) 
holdzl  =  int (grabedges [2, 0] +grabradius) 
holdz2  =  int (grabedges [2, 1] -grabradius) 
if  holdzl  ==  holdz2 : 

holdzl  =  holdzl  -  2 
holdz2  =  holdz2  +  2 
if  holdzl  >  holdz2: 

holdzl  =  int (grabedges [ 2 , 1 ] -grabradius ) 
holdz2  =  int (grabedges [ 2 , 0 ] fgrabradius ) 
x_c  =  float (random. randint (holdxl , holdx2 ) ) 
y_c  =  float (random. randint (holdyl , holdy2 ) ) 

z_c  =  float (random. randint (holdzl , holdz2 ) ) 

to 

dtree  =  sp  .  cKDTree  (refdat  [  :  ,  2  :  5]  ) 

neighpoint  =  dtree . query_ball_point ( [x_c, y_c, z_c] ,  grabradius) 
datsphere  =  ref dat [np . asarray (neighpoint ),  2 : 5 ] 
datspherel  =  datsphere [:, 0 : 3] - [x_c,y_c, z_c] 

R  =  np . array ([ [np . cos (angle )+ (axis [ 0 ] **2 )*( 1-np . cos (angle )), \ 

axis [0] *axis [1] *  ( 1-np .cos (angle) ) -axis [2] *np . sin (angle) , \ 
axis [0] *axis [2] *  ( 1-np .cos (angle) ) +axis [1] *np.sin (angle) ] , \ 
[axis [1] *axis [0] * (1-np .cos (angle) ) +axis [2] *np . sin (angle) , \ 
np .  cos (angle)  +  (axis [1] **2) * (1-np . cos  (angle) ) ,  \ 
axis [1] *axis [2] * (1-np. cos (angle) ) -axis [0] *np.sin (angle) ] , \ 
[axis [2] *axis [0] * (1-np .cos (angle) ) -axis [1] *np . sin (angle) , \ 
axis [2] *axis [1] *  (1-np. cos (angle) ) +axis [0] *np.sin (angle) , \ 
np .  cos (angle)  +  (axis [2] **2) *  (1-np . cos (angle) )  ] ] ) 
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datsphererot  =  np.dot (datspherel,  R) 
datspheref inal  =  datsphererot  +  grain_center 
return  datsphere,  datspheref inal 

def  check_ownership (grain_centers , atoms , box_limits ) : 

#checks  to  make  sure  an  atom's  nearest  center  is  its  parent  center,  deletes  atom  if  not 
grain_centers  =  ghost_centers (grain_centers, box_limits) 
centertree  =  sp . cKDTree (grain_centers ) 
dist,  label  =  centertree . query (atoms [:, 0 : 3] ) 
check  =  np . equal (atoms [:, 3 ], label ) 
output_atoms  =  atoms [check, 0 : 4 ] 
return  output_atoms ,  grain_centers 

to 

^  def  transpose_periodic_bounds (atoms , box_lims ) : 

#reflects  atoms  which  exist  outside  the  box  bounds  through  the  periodic  conditions 
lower  =  np. empty ( [len (atoms) , 3] ) 

lower[:,0:3]  =  [box_lims [ 0 , 0 ] , box_lims [ 1 , 0 ] , box_lims [2 , 0 ] ] 
upper  =  np. empty ( [len (atoms) , 3] ) 

upper[:,0:3]  =  [box_lims [ 0 , 1 ] , box_lims [ 1 , 1 ] , box_lims [2 , 1 ] ] 
test_l  =  np. less_equal (atoms [ : , 0 : 3] , lower [ : , 0 : 3] ) 
test_u  =  np . greater_equal (atoms [ : , 0 : 3 ]  ,  upper [ : ,  0 : 3 ] ) 
shift_l  =  np. where (test_l, 1, 0) 
shift_u  =  np . where (test_u, -1 , 0 ) 

shift_L  =  np .multiply ( shift_l [ : , 0 : 3 ] , np . array ( [box_lims [0,1] -box_lims [ 0 , 0 ] , \ 

box_l ims [1,1] -box_l ims [ 1 , 0 ] , \ 
box_l ims [2,1] -box_l ims [ 2 , 0 ] ] ) ) 


Approved  for  public  release;  distribution  is  unlimited. 


shift_U  =  np .multiply ( shift_u [ : , 0 : 3 ] , np . array ( [box_lims [0,1] -box_lims [ 0 , 0 ] , \ 

box_l ims [1,1] -box_l ims [ 1 , 0 ] , \ 
box_l ims [2,1] -box_l ims [ 2 , 0 ] ] ) ) 

out_atoms  =  atoms 

out_atoms [ : , 0 : 3 ]  =  np . add (atoms [ : , 0 : 3 ] , shift_L) 
out_atoms [ : , 0 : 3 ]  =  np . add (atoms [ : , 0 : 3 ] , shift_U) 
return  out_atoms 

def  ref erence_rescale (ref_atoms, box_lims, grain_size) : 

#function  resizes  the  reference  structure  to  accomodate  requested  grain  size 
rescale_length  =  np . array ([ (grain_size  -  (box_lims [ 0, 1 ] -box_lims [0, 0 ] ) ) , \ 

(grain_size  -  (box_lims [ 1 , 1 ] -box_lims [ 1, 0 ] ) ) , \ 

(grain_size  -  (box_lims [2, 1] -box_lims [2, 0] ) ) ] ) 

to 

^  if  rescale_length [ 0 ]  >0.0  \ 

or  rescale_length [ 1 ]  >  0.0\ 

or  rescale_length [2 ]  >  0.0: 

n_x  =  int (rescale_length [0] / (box_lims [0, 1] -box_lims [0, 0] ) ) +2 
n_y  =  int (rescale_length [ 1 ] / (box_lims [ 1 , 1 ] -box_lims [ 1 , 0 ] ) ) +2 
n_z  =  int (rescale_length [2 ] / (box_lims [2 , 1 ] -box_lims [ 2 , 0 ] ) ) +2 
new_ref  =  ref_atoms 
for  i  in  range ( 0 , n_x+l ) : 

for  j  in  range ( 0 , n_y+l ) : 

for  k  in  range ( 0 , n_z+l ) : 

if  i  ==  0  and  k  ==  0  and  j  ==  0: 
continue 


else : 
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K> 


else : 


new_ats  =  ref_atoms  +  [ 0 , 0 ,  float ( i )* (box_lims [ 0 , 1 ] -box_lims [ 0 , 0 ]), \ 

float ( j ) * (box_lims [1,1] -box_lims [ 1 , 0 ] ) , \ 
float (k) * (box_lims [2,1] -box_lims [2, 0] ) ] 
new_ref  =  np . concatenate ( (new_ref, new_ats ) ) 


replicate_range  =  np . array ( [ [box_lims [0,0] , box_lims [0,0] +2  *rescale_length [ 0 ] ] , \ 

[box_lims [1,0] , box_lims [1,0] +2  *rescale_length [ 1 ] ] , \ 
[box_lims [2,0] , box_lims [2,0] +2  *rescale_length [ 2 ] ] ] ) 
check_x=np . where (ref_atoms [ : , 2 ]  <=  replicate_range [ 0 , 1 ] , True, False) 
rep_atoms_x  =  ref_atoms [check_x [ : ]  ,  : ] 

rep_atoms_x [ : , 2 ] =rep_atoms_x [ : , 2 ]  +  (box_lims [0,1] -box_lims [0,0] ) 

rep_atoms  =  np . concatenate (( ref_atoms , rep_atoms_x) , axis=0 ) 
check_y=np . where ( rep_atoms [ : , 3 ]  <=  replicate_range [ 1 , 1 ] , True, False) 
rep_atoms_y  =  rep_atoms [check_y [ : ]  ,  : ] 

rep_atoms_y [ : , 3 ] =rep_atoms_y [ : , 3 ]  +  (box_lims [1,1] -box_lims [1,0] ) 
rep_atoms  =  np . concatenate (( rep_atoms , rep_atoms_y) , axis=0 ) 
check_z=np . where (rep_atoms [ : , 4 ]  <=  replicate_range [2 , 1 ] , True, False) 
rep_atoms_z  =  rep_atoms [check_z [ : ] , : ] 

rep_atoms_z [ : , 4 ] =rep_atoms_z [ : , 4 ]  +  (box_lims [2,1] -box_lims [2,0] ) 
rep_atoms  =  np . concatenate ( (rep_atoms, rep_atoms_z) , axis=0) 
new_ref  =  rep_atoms 
return  new_ref 


def  multireference (ref_list, grain_size, alimits, rescale, percent_twin, twin_reference, axis, angle) : 

#function  will  extract  centroid  and  box  data  from  config  file  and  build  a  new  structure  accordingly 
=  np . loadtxt ( ' %s '  %  (ref_list) , skiprows=2,  \ 


c_ref s 
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N> 

00 


dtype= { ' names ' : ( ' id ' , ' c_x ' , ' c_y ' , ' c_z ' , ' axisr_x ' , ' axisr_y ' , ' axisr_z ' , ' angler ' , ' file ' ) , \ 
' formats '  :  ( ' i4 ' ,  ' f 4 ' ,  ' f 4 ' ,  ' f 4 ' ,  ' f 4 ' ,  ' f 4 1 ,  ' f 4 ' ,  ' f 4 ' ,  ' S32  '  )  } ) 
n_grains  =  f loat ( len (c_ref s ) ) 
n_twinned  =  n_grains  *  percent_twin 
check  =  np . empty (( int (n_twinned) , 1 ), dtype=int ) 
print  'reassigning  %i  reference  files'  %  (n_twinned) 
for  i  in  range ( 0 , int (n_twinned) ) : 

itl  =  random. randint (0, len (c_refs) -1) 
while  np . any (check [ : ] ==itl ) ==True : 

itl  =  random . randint ( 0 , len (c_refs ) -1 ) 
c_ref s [ ' f ile ' ] [itl]  =  '%s'  %  (twin_reference) 
check  [i, 0]  =  itl 

center  =  np. empty ( (len (c_refs) , 3) ) 

f=open ( ' %s_centroids . conf ig '  %  (output_name) , ' w ' ) 
f . write (' #data  for  centroids\n ' ) 

f.write('%f  %f  %f  %f  %f  %f  %f\n'  %  (grain_size*rescale,  alimits [ 0 , 0 ] ,  alimits [ 0 , 1 ] *rescale, \ 
alimits [1, 0] ,  alimits [1, 1 ] *rescale,  alimits [2 , 0 ] ,  alimits [ 2 , 1 ] *rescale ) ) 
finaldat  =  np . empty ( [0 , 4 ] ) 
print  'populating  atoms' 
for  i  in  range (0, len (c_refs) ) : 

print  'populating  grain  %i '  %  (i) 

limsx  =  np . f romstring ( line cache . get line ( ' %s '  %  ( c_ref s ['file'] [ i ] ) , 6) , sep= '  ' ) 

limsy  =  np . fromstring (linecache . getline ( ' %s '  %  (c_ref s [ ' f ile ' ] [ i ] ) , 7 ) , sep= '  ') 

limsz  =  np . fromstring (linecache . getline ( ' %s '  %  ( c_ref s ['file'] [ i ] ) , 8 ) , sep= '  ' ) 

limits  =  np . concatenate ( (limsx, limsy, limsz ) , axis=0 )  . reshape  (3,2) 
cdat=np. loadtxt ( ' %s '  %  c_refs ['file'] [i], skiprows=9, usecols=(0,l,2,3,4) ) 
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K> 

so 


if  1 . 5*grain_size*rescale  >=  (limits [0, 1] -limits [0, 0] )  \ 

or  1 . 5*grain_size*rescale  >=  ( limits [ 1 , 1 ] -limit s [ 1 , 0 ] )  \ 

or  1 . 5*grain_size*rescale  >=  ( limits [ 2 , 1 ] -limit s [ 2 , 0 ]) : 
cdat  =  ref erence_rescale (cdat, limits, 1 . 5*grain_size*rescale) 

limits  =  np . array ( [ [np . amin (cdat [ : , 2 ] ) , np . amax (cdat [ : , 2 ] ) ] , [np . amin (cdat [ : , 3 ] ) , \ 
np . amax (cdat [  :  ,  3  ]  )  ]  , [np . amin (cdat [ : , 4 ] ) , np . amax (cdat [ : , 4 ] ) ] ] ) 
if  axis  ! =  ' none ' : 

c_ref s [ ' axisr_x ' ]  =  axis[0] 
c_ref s [ ' axisr_y ' ]  =  axis[l] 
c_ref s [ ' axisr_z ' ]  =  axis  [2] 
haxis  =  axis 

f.write('%i  %f  %f  %f  %f  %f  %f  %f  %s\n'  %  (i,  c_refs [ 'c_x' ] [i] *rescale,  \ 

c_ref s [ ' c_y ' ] [i] *rescale,  c_refs['c_z'] [i]*rescale,  c_refs [  axisr_x ' ] [i] , \ 
c_ref s [ ' axisr_y ' ] [ i ] ,  c_ref s [ ' axisr_z ' ] [ i ] ,  c_ref s [  angler ' ] [ i ] ,  \ 
c_r e  f  s [ ' f i 1 e ' ]  [i])) 

center [ i, : ] = [c_ref s [ ' c_x ' ] [ i ] * rescale, c_ref s [ ' c_y '] [  i  ]  * re scale, c_ref s [ ' c_z ' ] [ i ] * re scale] 
[trydat,  seeddat ] =seed_rotations ( cdat , [ c_ref s [ ' c_x ' ] [i] *rescale, \ 

c_ref s [ ' c_y ' ] [i] *rescale, \ 
c_refs [ ' c_z ' ] [i] *rescale] , \ 

1 . 5*rescale*grain_size/2 , limits, haxis, \ 
c_ref s [ ' angler ' ] [ i ] ) 

%  (i,  c_ref s [ ' c_x ' ] [ i ] *rescale,  \ 

c_ref s [ ' c_y ' ] [i]*rescale,  \ 
c_ref s [ ' c_z ' ] [i]*rescale,  \ 
c_ref s [ ' axisr_x ' ] [i],  \ 


else : 

f . write ( ' %i  %f  %f  %f  %f  %f  %f  %f  %s\n' 
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c_ref s [ ' axisr_y ' ] [i],  \ 

c_ref s [ ' axisr_z ' ] [i],  \ 

c_ref s [' angler ' ] [i],  \ 
c_ref s [ ' f ile  '  ]  [ i ]  )  ) 

center [ i, : ] = [c_ref s [ ' c_x ' ] [ i ] * rescale, c_ref s [ ' c_y '] [  i  ]  * re scale, c_ref s [ ' c_z ' ] [ i ] * re scale] 
haxis  =  np . array ( [ c_ref s [ ' axisr_x ' ] [ i ] ,  c_ref s [ ' axisr_y ' ] [ i ] ,  c_ref s [ ' axisr_z ' ] [ i ] ] ) 
[trydat,  seeddat ] =seed_rotations ( cdat , [ c_ref s [ ' c_x ' ] [i] *rescale, \ 

c_ref s [ ' c_y ' ] [i] *rescale, \ 
c_refs [ ' c_z ' ] [i] *rescale] , \ 

1 . 5*rescale*grain_size/2 , limits,  haxis,  \ 
c_ref s [ ' angler ' ]  [ i  ]  ) 


seeddat  =  np . insert ( seeddat , 3, i, axis=l ) 

O  finaldat  =  np . concatenate (( finaldat , seeddat ) ) 

linecache . clearcache ( ) 
f . close ( ) 

return  finaldat,  center 

# - Main  Body - 

if  ref erence_mode  ==  'single': 

cdat  =  np . loadtxt ( ' %s '  %  (crref) , skiprows=9,usecols= (0, 1, 2, 3, 4) ) 
limsx  =  np . fromstring ( linecache . get line (' %s '  %  (crref ) , 6) , sep= '  ') 
limsy  =  np . fromstring ( linecache . get line (' %s '  %  (crref ) , 7) , sep= '  ') 
limsz  =  np . fromstring ( linecache . get line (' %s '  %  (crref ) , 8) , sep= '  ') 
lims  =  np . concatenate (( limsx, limsy , limsz ), axis=0 ). reshape (3 , 2 ) 
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#  print  lims 

#  Check  if  the  average  grain  size  is  larger  than  the  box  bounds 

if  1 . 5*avg_grain_size  >=  (lims [0, 1] -lims [0, 0] )  or  1 . 5*avg_grain_size  >=  ( lims [ 1 , 1 ] -lims [ 1 , 0 ] )  \ 

or  1 . 5*avg_grain_size  >=  (lims [2 , 1 ] -lims [2 , 0 ] ) : 
cdat  =  ref erence_rescale (cdat, lims, 1 . 5*avg_grain_size) 
lims  =  np . array ([ [np . amin ( cdat [:, 2 ]), np . amax (cdat [:, 2 ])], \ 

[ np . amin (cdat [ : , 3 ] ) , np . amax ( cdat [ : , 3 ] ) ] , \ 

[ np . amin ( cdat [ : , 4 ] ) , np . amax ( cdat [ : , 4 ] ) ] ] ) 
lims2  =  np . f romstring (box_bounds , sep= ' , ' ) . reshape (3 , 2 ) 

if  avg_grain_size  >=  (lims2 [ 0, 1] -lims2 [ 0, 0] )  or  avg_grain_size  >=  ( lims2 [ 1 , 1 ] -lims2 [ 1 , 0 ] )  \ 

or  avg_grain_size  >=  (lims2 [2, 1] -lims2 [2, 0] ) : 
sys . exit (" grain  size  larger  than  defined  simulation  box,  change  box  sizes") 
print  'assigning  grain  centers' 
center  =  grain_centroid (lims2, avg_grain_size) 
finaldat  =  np . empty ( [ 0 , 4 ] ) 

a  =  np . ones (( len (cdat ), 1 )). reshape ( len (cdat ),  1 ) 
b  =  np . arange ( 1 , len ( cdat ) tl ). reshape (len ( cdat ),  1 ) 
cdat  =  np . concatenate ( (b, a, cdat [:, 2 : 5 ]), axis=l ) 
centout=np . empty ( (len (center) , 8) ) 

f=open ( ' %s_centroids . conf ig '  %  (output_name) , ' w ' ) 
f . write (' #data  for  centroids\n ' ) 

f.write('%f  %f  %f  %f  %f  %f  %f\n'  %  (avg_grain_size, lims2 [ 0 , 0 ] , lims2 [ 0 , 1 ] , lims2 [ 1 , 0 ] , \ 

lims2 [1,1], lims 2 [2,0], lims2 [2,1])) 

print  'populating  atoms' 

for  i  in  range (0, len (center) ) : 

print  'populating  grain  %i '  %  (i) 
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if  direction  ==  'none': 

axis  =  np . array ( [ random. uniform ( 0 , 10 ) , random . uniform ( 0 , 10 )  ,  random . uniform (0,10)]) 
axnorm  =  np . linalg . norm (axis ) 
axis  =  axis/axnorm 

angle  =  mt . radians ( random. uniform ( 0 , 180 ) ) 
else : 

axis  =  direction 

angle  =  random. uniform ( 0  ,  180 ) 

f.write('%i  %f  %f  %f  %f  %f  %f  %f  %s\n'  %  ( i, center [ i, 0 ], center [ i, 1 ], center [ i, 2 ], axis [ 0 ], \ 

axis [ 1 ] , axis [2 ] , angle, crref ) ) 

[trydat,  seeddat]  =  seed_rotations (cdat, center [i, 0 : 3] , 1 . 5*avg_grain_size/2, lims, axis, angle) 
seeddat  =  np . insert (seeddat, 3, i, axis  =  1) 
finaldat  =  np . concatenate (( finaldat , seeddat ) ) 
f . close ( ) 

print  'finaldat=',  len (finaldat) 
print  'trimming  excess  atoms' 

cut_atoms, check_centers  =  check_ownership (center, finaldat [ : , 0 : 4 ] , lims2) 
print  ' cut_atoms= ' ,  len (cut_atoms) 
print  'wrapping  periodic  bounds' 

fin_atoms  =  transpose_periodic_bounds (cut_atoms,  lims2) 
idcol  =  np . ones ([ len (cut_atoms ), 1 ] ) 

fin_atoms  =  np . concatenate ( (idcol, cut_atoms) ,  axis  =  1) 

numcol  =  np . arange ( 1 , len ( f in_atoms ) +1 ) . reshape (len ( f in_atoms ) , 1 ) 

fin_atoms  =  np . concatenate  (  (numcol, fin_atoms) ,  axis  =  1) 

output_atoms  =  np . column_stack ( ( f in_atoms [ : , 0 ] , f in_atoms [ : , 5 ] +1 , f in_atoms [ : , 2 ] , \ 
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f in_atoms [ : , 3 ] , f in_atoms [ : , 4 ] ) ) 


print  'outputting  to  data  file' 

head  =  ("#data  file  for  lammps\n%i  atoms\n%s  atom  types\n%f  %f  xlo  xhi\n%f  %f  ylo  "\ 

"yhi\n%f  %f  zlo  zhi\n\nAtoms\n" )  %  (len (output_atoms) , len (center) +1, lims2 [0, 0] , \ 

lims2 [0,1], lims2 [1,0], lims2 [1,1], lims2 [2,0], lims2 [2,1]) 
np . savetxt ( ' data . %s '  %  (output_name) , output_atoms,  fmt='%i  %i  %f  %f  %f') 
f=open ( ' data . %s '  %  (output_name) , ' r ' ) 
datal  =  f.readO 

f .  close ( ) 

g=open ( ' data . %s '  %  (output_name) , ' w ' ) 

g.  write (' %s\n '  %  (head)) 
g.write('%s'  %  (datal)) 
g. close () 

else : 

print  'reading  configuration  data' 

conf  =  np . f romstring ( linecache . getline ( ' %s '  %  (ref_list ) , 2 ) , sep= '  ') 

avg_grain_size  =  conf[0] 

lims  =  conf [1 :]. reshape (3, 2) 

lims2  =  lims  *  rescale_f actor 

axis  =  direction 

angle  =  mt . radians ( random . uniform ( 0 , 180 ) ) 

[finaldat,  center]  =  multireference (ref_list, avg_grain_size, lims, rescale_f actor , contwin, \ 

twinr, axis, angle) 

if  avg_grain_size*rescale_f actor  >=  (lims2 [0, 1] -lims2 [0, 0] )  \ 
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or  avg_grain_size*rescale_f actor  >=  ( lims2 [ 1 , 1 ] -lims2 [ 1 , 0 ] )  \ 

or  avg_grain_size*rescale_f actor  >=  (lims2 [2, 1] -lims2 [2, 0] ) : 
sys . exit (" grain  size  larger  than  defined  simulation  box,  change  box  sizes") 

a  =  np . ones (( len ( finaldat ), 1 )). reshape (len ( finaldat ), 1 ) 
b  =  np . arange ( 1 , len ( finaldat ) +1 ). reshape ( len ( finaldat ),  1 ) 
finaldat  =  np . concatenate ( (b, a, finaldat ), axis=l ) 

print  'trimming  excess  atoms' 

cut_atoms, check_centers  =  check_ownership (center, finaldat [ : , 2 : 6] , lims2) 
a  =  np . ones (( len (cut_atoms) , 1 )). reshape ( len (cut_atoms ), 1 ) 
b  =  np . arange ( 1 , len ( cut_atoms ) +1 ). reshape ( len (cut_atoms ), 1 ) 
cut_atoms  =  np . concatenate ( (b, a, cut_atoms ), axis=l ) 

print  'wrapping  periodic  bounds' 

fin_atoms  =  transpose_periodic_bounds (cut_atoms [ : , 2 : 5] , lims2) 
idcol  =  np . ones ([ len (cut_atoms ), 1 ] ) 

fin_atoms  =  np . concatenate ( (idcol, cut_atoms) ,  axis  =  1) 

numcol  =  np . arange ( 1 , len ( f in_atoms ) +1 ) . reshape (len ( f in_atoms ) , 1 ) 

fin_atoms  =  np . concatenate  (  (numcol, fin_atoms) ,  axis  =  1) 

output_atoms  =  np . column_stack ( ( f in_atoms [ : , 0 ] , f in_atoms [ : , 7 ] +1 , f in_atoms [ : , 4 ] , \ 

fin_atoms [ : , 5] , fin_atoms [ : , 6] ) ) 


print  'outputting  to  data  file' 

head  =  ("#data  file  for  lammps\n%i  atoms\n%s  atom  types\n%f  %f  xlo  xhi\n%f  %f  ylo  "\ 
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"yhi\n%f  %f  zlo  zhi\n\nAtoms\n" )  %  (len (fin_atoms) , len (center) +1,  lims2 [0,  0]  ,  \ 

lims2 [0,1], lims2 [1,0], lims2 [1,1], lims2 [2,0], lims2 [2,1]) 
np . savetxt ( ' data . %s '  %  (output_name) , output_atoms ,  fmt='%i  %i  %f  %f  %f') 
f=open ( ' data . %s '  %  (output_name) , ' r ' ) 
datal  =  f.readO 

f .  close ( ) 

g=open ( ' data . %s '  %  (output_name) , ' w' ) 

g.  write (' %s\n '  %  (head)) 
g.write('%s'  %  (datal)) 
g . close ( ) 
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